Efficient Nyström-type method for the solution of highly oscillatory Volterra integral equations of the second kind

Highly oscillatory Volterra integral equations are frequently encountered in engineering applications. The Nyström-type method is an important numerical approach for solving such problems. However, there remains scope to further optimize and accelerate the Nyström method. This paper presents a novel Nyström-type method to efficiently approximate solutions to second-kind Volterra integral equations with highly oscillatory kernels. First, the unknown function is interpolated at Chebyshev points. Then the integral equation is solved using the Nyström-type method, which leads to a problem of solving a system of linear equations. A key contribution is the technique to express the fundamental Lagrange polynomial in matrix form. The elements of the matrix, which involves highly oscillatory integrals, are calculated by using the classical Fejér quadrature formula with a dilation technique. The proposed method is more efficient than the one proposed in the recent literature. Numerical examples verify the efficiency and accuracy of the proposed method.


Introduction
Highly oscillatory Volterra integral equations of the second kind, which can be described as follows, are frequently encountered in several engineering applications, including wave propagation, quantum mechanics, acoustic scattering, signal processing, and image recognition.
where y(x) is the unknown function, f(x) is a given function defined on [a, b], and K(x, t) is a known highly oscillatory function defined on D = {(x, t): a � t � x � b}, which depends on a real positive fixed parameter ω.Examples of such kernels are Kðoðx À tÞÞ ¼ expðioðx À tÞÞ and Kðoðx À tÞÞ ¼ J n ðoðx À tÞÞ; where J ν is the Bessel function of the first kind of order ν � 0. When w � 1, the kernel function is highly oscillatory.Discretizing the integral equation using the collocation or Galerkin methods will involve the calculation of the integral of a highly oscillatory function.Collocation and Galerkin methods based on traditional quadrature formulae are computationally costly, which makes it extremely difficult to solve the integral Eq (1) numerically.The computation of the integrals of kernel functions and unknown functions, along with the discretization of integral equations, are the main challenges in solving integral equations.To address this problem, many methods have been proposed for the calculation of integrals of highly oscillatory functions.These methods include asymptotic methods, Filon-type method, Levin-type method, numerical steepest descent and generalized quadrature rules [1][2][3][4][5][6][7][8].All these methods have the feature that the higher the frequency, the higher the calculation accuracy.But for general oscillators, the above methods are sometimes complicated.
In order to meet the demand of the application of highly oscillatory integral equations, many scholars have proposed some numerical methods for solving these kinds of integral equations by combining the methods for calculating the highly oscillatory integrals with the collocation methods and the discontinuous Galerkin et al., and providing the error estimates of the numerical methods.In 2004, Davies and Duncan derived the stability and convergence of the solution of the collocation method for the first kind of Volterra integral equation with highly oscillatory Bessel kernel [9].In 2009, Brunner et al. investigated the discontinuous Galerkin method for solving a class of highly oscillatory integral equations.In 2011, Wang and Xiang et al. studied the numerical solution of a class of highly oscillatory Volterra integral equations by combining asymptotic methods and Filon-type methods for the calculation of highly oscillatory integrals, and gave an error analysis [10].In 2013, Xiang and Brunner studied the collocation method to solve the weakly singular highly oscillatory integral equation with Bessel kernel, analyzed some asymptotic properties of the solution of the equation, and proposed an efficient algorithm to solve the integral equation numerically [11].Meanwhile, Xiang and Wu investigated the numerical solution of a class of Volterra integral equations with highly oscillatory triangular kernel by combining the Filon method and the collocation method, and gave an error analysis [12].In the same period, Xiang and He investigated the discontinuous Galerkin method for solving this class of integral equations and found that the computational accuracy of the discontinuous Galerkin method is comparable to that of the collocation method, but the collocation method is simpler [13].In 2014, Wu studied the case of solving the weakly singular highly oscillatory Volterra integral equation by the collocation method on graded meshes.It was found that for weak singularities, the graded meshes is effective at low frequencies, while in the case where the kernel function is highly oscillatory, the employed graded meshes are not very useful [14].Meanwhile, Brunner studied the properties of highly oscillatory Volterra integral operators from a theoretical point of view [15].In 2015, Ma and Fang et al. studied the direct Filon method for the same kind of integral equations and gave a better error estimation formula [16].In 2015, Brunner H., Ma Y., Xu Y. studied the highly oscillatory nature of the solutions of highly oscillatory integral differential equations [17].In 2018, Xiang S et al. proposed a numerical solution to a class of highly oscillatory delay potential integral equations by combining the collocation method with the traditional fast multipole method [18].In 2019, Zhang Q., Xiang S. combined the collocation method with the traditional fast multipole method to present the numerical solution of a class of highly oscillatory Volterra integral equations and gave error estimates.And it concluded that the fast multipole method can greatly improve the efficiency of the collocation method as the number of collocation points increases [19].In 2021, using the Block-FMM method based on Chebyshev polynomial expansion, Wu et al. proposed an efficient method for solving the Volterra integral equation and gave error estimates [20].In 2022, Longbin Zhao et al. derived a method for solving the highly oscillatory Volterra integral equation and analyzed the asymptotic behavior of the solution to obtain an error estimation formula with respect to both frequency and step size [21,22].
Recently, by combining the Legendre polynomials and Gauss-Legendre quadrature formulas, Luisa et al. presented a Nystro ¨m-type method for the numerical approximation of the solution of Volterra integral Eq (1), and the convergence of the method was derived.This new numerical method does not use the method of computing highly oscillatory integrals [23].
In this paper, we combine Chebyshev polynomials with the Nystro ¨m-type method to study the numerical solution of the following Volterra integral equations of second kind with highly oscillatory kernel where y(x) is the unknown function, f(x) and K(ω(x − t)) are given functions.For the general interval [a, b], it can be mapped to the interval [−1, 1] by the transformation

The advantages of the method
• Representing Lagrange polynomials in matrix form improves computational efficiency and simplifies coding.
• The Feje ´r quadrature with dilation avoids specialized oscillatory integral techniques while maintaining accuracy.
• The method maintains high accuracy using standard numerical integration.
• The method enhances computational speed and ease of implementation.

The disadvantages of the method
• Increasing the parameter ω requires more quadrature points, reducing efficiency.
• Extremely oscillatory solutions decrease the accuracy due to limited quadrature resolution.
In summary, the matrix Lagrange basis and Feje ´r quadrature offer advantages in accuracy, speed, and coding simplicity.The main limitations occur for very large ω values, where efficiency degrades, and for solutions with sharp oscillations, where accuracy declines.
This paper is organized as follows: In Section 2, we introduce the basic concepts and formulas of Chebyshev polynomials and the Fej'er quadrature formula.In Section 3, we derive a method to solve integral Eq (2).In Section 4, we derive convergence rates for the proposed method.In Section 5, we present some numerical examples to demonstrate the efficiency and accuracy of the proposed method.

Chebyshev polynomials and Feje ´r quadrature formula
In this section, the basic properties of Chebyshev polynomials and the Feje ´r quadrature formula are introduced, for details see [24,25].

Chebyshev polynomials and the Lagrange interpolation process
The Chebyshev polynomial T n (x) of the first kind of degree n is defined by the relation T n ðxÞ ¼ cosðnyÞ; when x ¼ cosðyÞ: Using the trigonometric identity the fundamental recurrence relation is obtained T n ðxÞ ¼ 2xT nÀ 1 ðxÞ À T nÀ 2 ðxÞ; n ¼ 2; 3; . . .: The initial conditions efficiently generate all T n (x) recursively.
The Chebyshev polynomial T n (x) has n roots located at which are commonly referred to as the Chebyshev nodes.Denote fx k g n k¼1 as the zeros of T n (x), the approximating Lagrange polynomial based on these nodes can be expressed as where ℓ k is the k-th fundamental Lagrange polynomial defined as On the other hand, the approximating polynomial L n ð f ; xÞ to the function f(x) can be expressed as where : and x m are the roots of T n (x).By changing the order of summation, L n ð f ; xÞ can be written in the form of where Using the matrix form can improve computational efficiency.

Feje ´r quadrature formula
This subsection introduces basic concepts and formulas for interpolatory integration formulas using Chebyshev nodes.For more details, see [25].For a continuous function f on [−1, 1], we can approximate the integral by sums.Let This integral can be approximated as For interpolatory rules using Chebyshev nodes, w k is given by Using the Christoffel-Darboux formula for T n (x) gives Eq (3) can be written as From T 0 n ðcos yÞ ¼ nðsin nyÞ=sin y, we have ( Substituting these formulas into (4) leads to Therefore, Feje ´r's integration formula is where w k defined by (5).

Nystro ¨m-type method for Volterra integral equations
Define operator V : C 0 ð½À 1; 1�Þ !C 0 ð½À 1; 1�Þ as then the Eq (2) can be written in the operator form To represent the unknown function y(x), we define the Lagrange polynomial based on Chebyshev nodes as follows By substituting this expression into (7), we obtain the following equation When collocating at the zeros x k of the Chebyshev polynomials, we derive the linear system of equations Here, y j = y(x j ) represents the unknowns, and the coefficients c j (x k ) are defined as follows S n ðx j ; tÞKðoðx À tÞÞdt; j ¼ 1; 2; . . .; n: ð8Þ The preceding system system of order n × n can be expressed in matrix form as follows , and f ¼ ½f ðx k Þ� n k¼1 .Once the linear system of Eq (9) has been solved, we can calculate the Nystro ¨m interpolant as follows In general, the coefficients c j (x k ) cannot be computed analytically.The next question is how to efficiently calculate c j (x k ).In [23], a method known as the "dilated product quadrature formula" is introduced.This method involves expanding the interval and subdividing it into numerous smaller intervals to ensure that the integrand function on each small interval is non-oscillatory.We apply the same technique and consider Feje ´r quadrature formula.
With the change of variable t = z/ω, we can express Dividing the interval into N subintervals each of length δ gives By remapping each integral to [−1, 1] with the transformation each integral can be approximated using the n-point Feje ´r quadrature formula.

Convergence and error analysis
By utilizing the asymptotic behavior of the coefficients in the Chebyshev series expansion of function and by considering the aliasing errors of Feje ´r's first and second rules for approximating integrals, Xiang derive improved convergence rates for Feje ´r's rules.The key innovations involve utilizing asymptotics to derive tighter bounds on aliasing errors and combining this information with existing error analysis to achieve overall sharper rates, as stated in the following lemma.
Lemma 1 [26] Suppose f(x) has an absolutely continuous where I½ f � ¼ R 1 À 1 f ðxÞdx, I n [f] denotes the (n + 1)-point Fejér's first rules and zðk þ 1Þ ¼ P 1 p¼1 1 p kþ1 denotes the Riemann zeta function.Let c j (x) and c n j ðxÞ be the coefficients defined in ( 8) and (11), respectively.Assume that S n x; φ À 1 i ðtÞ o � � Kðox À φ À 1 i ðtÞÞ has an absolutely continuous (k − 1)st derivative on [−1, 1] and a kth derivative of bounded variation V k for some k � 1.Then, Define ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ð1 À tÞðt þ 1Þ p : The corresponding approximate equation can then be written as Denote je n m ðy; xÞj ¼ jðVyÞðxÞ À ðV n m yÞðxÞj � jðVyÞðxÞ À ðV m yÞðxÞj þ jðV m yÞðxÞ À ðV n m yÞðxÞj: Assuming that the kernel function K satisfies the following conditions jKðx; tÞj ffi ffi ffi ffi ffi ffi ffi ffiffi �ðtÞ Also, By applying Theorem 5.1.12in [27] we can obtain where E m ðyÞ ¼ inf P m 2P m ky À P m k 1 and P m denotes the space of the algebraic polynomials of degree at most m.Theorem 1 Denote by y the unique solution of Eq (7) and by y n m the unique solution of (14).Then, if f(x) 2 C r ([−1, 1]) and the kernel K satisfies the conditions (15) and (16), the following error estimate holds true where the constants in O do not depend on n and m.Proof From ( 7) and ( 14) we deduce that so that by applying (17) and ( 13) and the assumptions we obtain the result.

Numerical examples
This section provides several numerical examples to demonstrate the effectiveness and accuracy of the proposed algorithm.For each example, c j (x) is computed using the 20-point Feje ´r quadrature formula, and the solution is approximated using the Nystro ¨m interpolant in (10).The absolute error ε 20  n ðxÞ ¼ jyðxÞ À ŷ20 n ðxÞj is used to measure the algorithm's accuracy, where y(x) is the exact solution.For Example 2 and 3, we consider the precise approximations generated by the product rule with m = 750 and calculate the errors ε 20 750;n ðsÞ ¼ jy 20 750 ðxÞ À y 20 n ðxÞj: Example 1.Consider the following integral equation where y(x) = e x is the exact solution and the right-hand side f(x) is given by Since the kernel function, right-hand side, and solution are all analytic and non-oscillatory for Example 1, convergence is very rapid as seen in Table 1 and Fig 1 .For the same n and ω values, the computing time (e.g.0.1549s for n = 16, ω = 10 4 ) is approximately 10 times less than the 1.5381s reported in [23], demonstrating the improved efficiency of the proposed method.The calculated condition number of the matrix I n þ A n also indicates the system is extremely well-conditioned.
Example 2. Let us consider the following integral equation whose exact solution is not known.Example 3. Let us consider the following integral equation whose exact solution is not known.

Conclusion
This paper proposes a numerical method for solving highly oscillatory Volterra integral equations.The method combines Chebyshev polynomials, the Nystro ¨m approach, and n-point Feje ´r quadrature with dilation, avoiding reliance on specialized techniques for oscillatory integrals.A key efficiency enhancement is the representation of fundamental Lagrange polynomials in matrix form, improving computational speed and simplifying coding implementation.Numerical results demonstrate that the method retains high accuracy using standard numerical integration, despite the oscillatory nature of the problems.Experiments verify the effectiveness of the proposed technique on test examples.Compared to recent literature, the method provides improved computational efficiency.In summary, this work presents an accurate and computationally advantageous scheme for oscillatory Volterra integral equations.By leveraging these numerical techniques, the proposed Nystro ¨m-type method is efficient, straightforward to code, and suitable for practical applications.

For example 2 ,
Fig 2 illustrates that the kernel function exhibits increased oscillations as ω increases.Tables 2 to 3 and Figs 3 to 4 demonstrate the convergence of the method.It is easy to see that the scheme is efficient.Moreover, the condition number of the system does not depend on n.